#include "flow.h"
double rho_ref, u_ref, p_ref, vol_ref;
size_t flow::imax{0};
size_t flow::ib2{0};
/**
 * @brief 开辟内存空间
 * 
 * @param imax 包含 ghost 单元在内的单元数
 */
void flow::iniMemory(size_t imax) noexcept
{
    flow::imax = imax;
    flow::ib2 = imax - 1;

    cv.reshape(imax, ncVar);
    cvold.reshape(imax, ncVar);
    p.resize(imax);
    dt.resize(imax);
    diss.reshape(imax, ncVar);
    rhs.reshape(imax, ncVar);
    ls.reshape(imax, ncVar);
    rs.reshape(imax, ncVar);

    buffer1d.resize(imax);
    // cusp 格式有些特殊要求
    buffer2d.reshape(imax + 1, ncVar);
}

/**
 * @brief 初始化流场
 *      使用完全气体方程， 条件：入口总压、总温，出口静压
 *      将管内静压全部为出口静压，再利用总温和总压计算
 * @param area 
 */
void flow::iniFlow(const vec1d<double> &area) noexcept
{
    using inp::cpgas;
    using inp::gam1;
    using inp::gamma;
    using inp::gap1;
    using inp::p01;
    using inp::p2;
    using inp::rgas;
    using inp::t01;

    // iniMemory()
    // double gam1{gamma - 1.0}, gap1{gamma + 1.0};
    // double rgas{cpgas * gam1 / gamma};               ///气体常量 R
    double temp = t01 * pow(p2 / p01, gam1 / gamma); ///< 静温
    double rho{p2 / (rgas * temp)};
    double mach = sqrt(2 * (t01 / temp - 1.0) / gam1);
    double as = sqrt(gamma * p2 / rho); ///< 音速
    double u{as * mach};
    double massFlow{rho * u * area(1)};
    double E{(cpgas - rgas) * t01};

    size_t imax = area.size();
    for (size_t i = 0; i < imax; ++i)
    {
        cv(i, 0) = rho * area(i);
        cv(i, 1) = massFlow;
        cv(i, 2) = rho * E * area(i);
        p(i) = p2;
        // std::cout << cv(i,1) << "\n";
    }
    // exit(EXIT_FAILURE);
    // 限制器参考量？？？ 需要无量纲量
    rho_ref = rho;
    u_ref = u;
    p_ref = p2;
    vol_ref = 1.0;
}
/**
 * @brief 续算读入之前的解
 * 
 */
void flow::rsolut() noexcept
{
}

/**
 * @brief 边界条件
 *      入口亚音速、出口亚音速或超音速： 使用黎曼不变量设置
 * 
 * @param area 单元中心面积
 * @param ib2    
 */
void flow::bcond(const vec1d<double> &area, const size_t ib2) noexcept
{
    using inp::cpgas;
    using inp::gam1;
    using inp::gamma;
    using inp::gap1;
    using inp::p01;
    using inp::p2;
    using inp::rgas;
    using inp::t01;
    using std::cerr;

    // inlet (subsoinc)  268 269 页

    double ud{cv(1, 1) / cv(1, 0)};
    double cs2{gamma * p(1) / cv(1, 0) * area(1)}; ///< 音速平方
    double c02{cs2 + 0.5 * gam1 * ud * ud};
    // outgoing riemann invariant
    double rinv = ud - 2.0 * sqrt(cs2) / gam1;
    double dis = gap1 * c02 / (gam1 * rinv * rinv) - 0.5 * gam1;
    if (dis < 0.0)
    {
        cerr << "warning(bcond0): discriminant < 0 at inflow boundary \t"
             << dis << "\n";
        dis = 1.0e-20;
    }

    double cb = -rinv * gam1 / gap1 * (1.0 + sqrt(dis));
    double cc02 = (cb * cb / c02 < 1.0) ? cb * cb / c02 : 1.0;
    double tb{t01 * cc02};
    double pb = p01 * pow(tb / t01, gamma / gam1);
    double rhob{pb / (rgas * tb)};
    double ub = sqrt(2.0 * cpgas * (t01 - tb));
    // std::cout << std::setprecision(10) << ud << "\t" << cs2 << "\t" << c02 << "\t"
    //         << rinv << "\t" << dis << "\n";
    // std::cout << cb << "\t" << cc02 << "\t" << tb << "\t"
    //         << pb << "\t" << rhob << "\t" << ub << "\n";

    cv(0, 0) = rhob * area(1);
    cv(0, 1) = rhob * area(1) * ub;
    cv(0, 2) = area(1) * (pb / gam1 + 0.5 * rhob * ub * ub);
    p(0) = pb;

    // std::cout << cv(0, 0) << "\t" << cv(0, 1) << "\t"
    //         << cv(0, 2) << "\t" << p(0) << "\n";
    // outlet  subsonic and supersonic

    double rho{cv(ib2 - 1, 0) / area(ib2 - 1)};
    double u{cv(ib2 - 1, 1) / cv(ib2 - 1, 0)};
    double cs = sqrt(gamma * p(ib2 - 1) / rho);

    if (u < cs)
    { // 亚音速流
        pb = p2;
        rhob = rho + (pb - p(ib2 - 1)) / (cs * cs);
        ub = u + (p(ib2 - 1) - pb) / (rho * cs);
    }
    else
    { // 超音速流
        pb = p(ib2 - 1);
        rhob = rho;
        ub = u;
    }
    // imax = ib2+1
    cv(ib2, 0) = rhob * area(ib2);
    cv(ib2, 1) = rhob * area(ib2) * ub;
    cv(ib2, 2) = area(ib2) * (pb / gam1 + 0.5 * rhob * ub * ub);
    p(ib2) = pb;

    // std::cout << cv(ib2, 0) << "\t" << cv(ib2, 1) << "\t"
    //         << cv(ib2, 2) << "\t" << p(ib2) << "\n";
    // exit(EXIT_FAILURE);
}

/**
 * @brief 计算时间步长
 *      谱半径计算有些出人意料
 * 
 * @param vol 
 * @param area 
 * @param x 
 */
void flow::timeStep(const vec1d<double> &vol, const vec1d<double> &area,
                    const vec1d<double> &x) noexcept
{
    using inp::gamma;
    double cs, dx, sprad;
    for (size_t i = 1; i < ib2; ++i)
    {
        cs = sqrt(gamma * p(i) / cv(i, 0) * area(i));
        dx = 0.5 * (x(i + 1) - x(i - 1));
        sprad = cs * sqrt(dx * dx + area(i) * area(i)) + area(i) * std::abs(cv(i, 1) / cv(i, 0));
        dt(i) = vol(i) / sprad;
        // std::cout << std::setprecision(16) << std::scientific
        //         << dt(i) << "\t" << vol(i) << "\n";
        //         << diss(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
    dt(0) = dt(1);
    dt(imax - 1) = dt(imax - 2);
}

/**
 * @brief 计算人工粘性项
 *      注意小小的外推和求值方式
 * @param vol 
 * @param beta 
 */
void flow::dissip(const vec1d<double> &vol, double beta) noexcept
{
    using inp::vis2;
    using inp::vis4;
    vec1d<double> &p_sensor = buffer1d;
    vec2d<double> &prediss = buffer2d;
    // vec2d<double>
    for (size_t i = 1; i < ib2; ++i)
    {
        p_sensor(i) = std::abs(p(i + 1) - 2 * p(i) + p(i - 1)) / (p(i + 1) + 2 * p(i) + p(i - 1));
        // std::cout << std::setprecision(16) << std::scientific << p_sensor(i) << "\t"
        // << std::abs(3.5) << "\n";
    }
    // exit(EXIT_FAILURE);
    p_sensor(0) = p_sensor(1);
    p_sensor(imax - 1) = p_sensor(imax - 1);

    // disspation fluxes at i+1/2
    size_t im1, ip1, ip2;
    double eps2, eps4, eval, p_sensor_max;
    // 注意， unsinged 发生了溢出
    for (int i = 0; i < ib2; ++i)
    {
        // std::cout << i << "\n";
        // max
        im1 = (i - 1 < 0) ? 0 : i - 1;
        ip1 = i + 1;
        // min
        ip2 = (i + 2 < (imax - 1)) ? i + 2 : (imax - 1);

        eval = 0.5 * (vol(i) / dt(i) + vol(ip1) / dt(ip1));

        p_sensor_max = (p_sensor(i) < p_sensor(ip1)) ? p_sensor(ip1) : p_sensor(i);

        eps2 = vis2 * p_sensor_max * eval;
        eps4 = vis4 * eval;
        eps4 = (0.0 < eps4 - eps2) ? eps4 - eps2 : 0.0;

        prediss(i, 0) = eps2 * (cv(ip1, 0) - cv(i, 0)) -
                        eps4 * (cv(ip2, 0) - 3 * cv(ip1, 0) +
                                3 * cv(i, 0) - cv(im1, 0));

        prediss(i, 1) = eps2 * (cv(ip1, 1) - cv(i, 1)) -
                        eps4 * (cv(ip2, 1) - 3 * cv(ip1, 1) +
                                3 * cv(i, 1) - cv(im1, 1));

        prediss(i, 2) = eps2 * (cv(ip1, 2) - cv(i, 2)) -
                        eps4 * (cv(ip2, 2) - 3 * cv(ip1, 2) +
                                3 * cv(i, 2) - cv(im1, 2));
    }

    // hybrid multistage schemes
    double beta1{1.0 - beta};
    for (size_t i = 1; i < ib2; ++i)
    {
        diss(i, 0) = beta * (prediss(i, 0) - prediss(i - 1, 0)) + beta1 * diss(i, 0);
        diss(i, 1) = beta * (prediss(i, 1) - prediss(i - 1, 1)) + beta1 * diss(i, 1);
        diss(i, 2) = beta * (prediss(i, 2) - prediss(i - 1, 2)) + beta1 * diss(i, 2);

        // std::cout << std::setprecision(16) << std::scientific
        //         << diss(i, 0) << "\t" << diss(i,1) << "\t"
        //         << diss(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
}

/**
 * @brief 计算对流通量和残差
 *      守恒变量的平均
 * 
 * @param area 
 */
void flow::fluxCen(const vec1d<double> &area) noexcept
{
    // 守恒变量的平均（还有压强）， 其它一切物理量均由其推导
    double rhoav, rhouav, rhoEav;
    double rhoHav, pav, areaav, qs;
    vec2d<double> &flux = buffer2d;

    for (size_t i = 0; i < ib2; ++i)
    {
        areaav = 0.5 * (area(i) + area(i + 1));

        rhoav = 0.5 * (cv(i, 0) / area(i) + cv(i + 1, 0) / area(i + 1));
        rhouav = 0.5 * (cv(i, 1) / area(i) + cv(i + 1, 1) / area(i + 1));
        rhoEav = 0.5 * (cv(i, 2) / area(i) + cv(i + 1, 2) / area(i + 1));

        pav = 0.5 * (p(i) + p(i + 1));
        rhoHav = rhoEav + pav;
        qs = rhouav * areaav / rhoav;

        flux(i, 0) = rhoav * qs;
        flux(i, 1) = rhouav * qs + pav * areaav;
        flux(i, 2) = rhoHav * qs;
    }

    // flux + dissipation = rhs
    for (size_t i = 1; i < ib2; ++i)
    {
        rhs(i, 0) = flux(i, 0) - flux(i - 1, 0) - diss(i, 0);
        rhs(i, 1) = flux(i, 1) - flux(i - 1, 1) - diss(i, 1);
        rhs(i, 2) = flux(i, 2) - flux(i - 1, 2) - diss(i, 2);
        // std::cout << std::setprecision(16) << std::scientific
        //         << rhs(i, 0) << "\t" << rhs(i,1) << "\t"
        //         << rhs(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
}

void flow::SrcTerm(const vec1d<double> &area) noexcept
{
    double da;
    for (size_t i = 1; i < ib2; ++i)
    {
        da = 0.5 * (area(i + 1) - area(i - 1));
        rhs(i, 1) -= p(i) * da;
        // std::cout << std::setprecision(16) << std::scientific
        //         << rhs(i, 1) << "\n";
        //         << diss(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
}

/**
 * @brief 残值光顺
 * 
 */
void flow::Irsmoo() noexcept
{
    using inp::epsirs;

    vec1d<double> &d = buffer1d;
    d(0) = 0.0;
    rhs(0, 0) = 0.0;
    rhs(0, 1) = 0.0;
    rhs(0, 2) = 0.0;

    double t, eps2{2.0 * epsirs + 1.0};
    for (size_t i = 1; i < ib2; ++i)
    {
        t = 1.0 / (eps2 - epsirs * d(i - 1));
        d(i) = t * epsirs;
        rhs(i, 0) = t * (rhs(i, 0) + epsirs * rhs(i - 1, 0));
        rhs(i, 1) = t * (rhs(i, 1) + epsirs * rhs(i - 1, 1));
        rhs(i, 2) = t * (rhs(i, 2) + epsirs * rhs(i - 1, 2));
        // std::cout << std::setprecision(16) << d(i) << "\n";
        // std::cout << std::setprecision(16) << rhs(i, 0) << "\n";
    }

    size_t i = ib2 - 1;
    for (size_t j = 3; j <= ib2; ++j)
    {
        --i;
        rhs(i, 0) = rhs(i, 0) + d(i) * rhs(i + 1, 0);
        rhs(i, 1) = rhs(i, 1) + d(i) * rhs(i + 1, 1);
        rhs(i, 2) = rhs(i, 2) + d(i) * rhs(i + 1, 2);
        // std::cout << std::setprecision(16) << rhs(i, 0) << "\n";
    }
}

void flow::UpdateConservativeVar(const vec1d<double> &area) noexcept
{
    using inp::gam1;

    for (size_t i = 1; i < ib2; ++i)
    {
        cv(i, 0) = cvold(i, 0) - rhs(i, 0);
        cv(i, 1) = cvold(i, 1) - rhs(i, 1);
        cv(i, 2) = cvold(i, 2) - rhs(i, 2);

        p(i) = gam1 * (cv(i, 2) / area(i) - 0.5 * cv(i, 0) / area(i) * (cv(i, 1) / cv(i, 0) * cv(i, 1) / cv(i, 0)));
    }
}

void flow::Conver(double &drho, const string &convergence_file_name, const vec1d<double> &area) noexcept
{
    using inp::gamma;

    drho = 0.0;
    double drmax{0.0}, av_mass_flow{0.0};
    double dr{0.0}, rho{0.0}, u{0.0}, c{0.0};
    size_t index_drmax{0}, nsup{0};

    for (size_t i = 1; i < ib2; ++i)
    {
        dr = cv(i, 0) - cvold(i, 0);
        av_mass_flow += cv(i, 1);
        drho += dr * dr;
        if (std::abs(dr) > drmax)
        {
            drmax = std::abs(dr);
            index_drmax = i;
        }

        rho = cv(i, 0) / area(i);
        u = cv(i, 1) / cv(i, 0);
        c = sqrt(gamma * p(i) / rho);
        if (c < u)
            ++nsup;
    }
    av_mass_flow /= (ib2 - 1);

    static double drho1{0.0};
    static size_t iter{0};
    ++iter;
    // 以第一步平均密度差值作为基准
    if (iter == 1)
        drho1 = sqrt(drho / (ib2 - 1));

    drho = sqrt(drho / (ib2 - 1)) / drho1;

    std::cout << std::scientific << std::setprecision(9)
              << iter << "\t" << log10(drho) << "\t"
              << drmax << "\t" << index_drmax << "\t"
              << av_mass_flow << "\t" << cv(1, 1) - cv(ib2 - 1, 1) << "\t"
              << nsup << "\n";
    // 不应该这样干
    // using inp::convergence_file;
    // std::ofstream convergence_file{inp::fnconv, std::ios_base::out};
    // convergence_file << std::scientific << std::setprecision(9)
    //                  << iter << "\t" << log10(drho) << "\t"
    //                  << drmax << "\t" << index_drmax << "\t"
    //                  << av_mass_flow << "\t" << cv(1, 1) - cv(ib2 - 1, 1) << "\t"
    //                  << nsup << "\n";
    // // 不应该这样做， 先这样
    // convergence_file.close();
}

void flow::solver(const Grid &grid, double &drho) noexcept
{
    // using namespace inp;
    using inp::ark;
    using inp::betrk;
    using inp::cfl;
    using inp::epsirs;
    using inp::kdissip;
    using inp::ldiss;
    using inp::lsmoo;
    using inp::nrk;

    // store previous solution, set dissipation = 0;
    cvold.copyFrom(cv);
    diss = 0.0;
    // calculate time step
    timeStep(grid.vol, grid.area, grid.x);
    // std::cout << "4\n";
    double fac;
    const vec1d<double> &vol = grid.vol;

    for (size_t i = 0; i != nrk; ++i)
    {
        // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
        if (kdissip == "c" || kdissip == "C")
        {
            // do something
            // std::cout << "1\n";

            if (ldiss[i])
                dissip(vol, betrk[i]);
            // std::cout << "6\n";
            // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
            fluxCen(grid.area);
            // std::cout << "5\n";
            // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
        }
        else if (kdissip == "r" || kdissip == "R")
        {
            LRStateRoe(grid.area, grid.vol);
            FluxRoe(grid.area);
        }
        else
        {
            std::cerr << "now, only support central scheme and roe scheme\n";
            exit(EXIT_FAILURE);
        }

        SrcTerm(grid.area);
        // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";

        // 每一阶 残差前的系数
        fac = ark[i] * cfl;
        double adtv;
        // std::cout << fac << "\n";
        for (size_t i = 1; i < ib2; ++i)
        {
            adtv = fac * (dt(i) / vol(i));
            rhs(i, 0) *= adtv;
            rhs(i, 1) *= adtv;
            rhs(i, 2) *= adtv;
            // std::cout << std::setprecision(15) << std::scientific << i << "\t"
            //           << rhs(i, 0) << "\n";
        }
        // std::exit(EXIT_FAILURE);
        // std::cout << "sss\n";

        // 残差光顺
        if (lsmoo[i] && (epsirs > 0.0))
            Irsmoo();
        // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
        // exit(EXIT_FAILURE);
        // update
        UpdateConservativeVar(grid.area);

        bcond(grid.area, ib2);
    }

    Conver(drho, inp::fnconv, grid.area);
}

void flow::LRStateCusp(const vec1d<double> &area) noexcept
{
    vec2d<double> &du = buffer2d;
    // 注意， 移动了一位
    for (int i = 1; i < imax; ++i)
    {
        du(i, 0) = cv(i, 0) / area(i) - cv(i - 1, 0) / area(i - 1);
        du(i, 1) = cv(i, 1) / cv(i, 0) - cv(i - 1, 1) / cv(i - 1, 0);
        du(i, 2) = p(i) - p(i - 1);
    }
    du(0, 0) = du(1, 0);
    du(0, 1) = du(1, 1);
    du(0, 2) = du(1, 2);
    // buffer2d imax+1, 3
    du(imax, 0) = du(ib2, 0);
    du(imax, 1) = du(ib2, 1);
    du(imax, 2) = du(ib2, 2);

    auto CuspLimiter = [](double d1, double d2) {
        return 1.0 - (std::pow((d1 - d2) / (std::abs(d1) + std::abs(d2) + 1.0e-20), 2.0));
    };

    for (int i = 0; i < imax; ++i)
    {
        ls(i, 0) = cv(i, 0) / area(i) + 0.25 * (du(i, 0) + du(i + 2, 0)) * CuspLimiter(du(i, 0), du(i + 2, 0));
        ls(i, 1) = cv(i, 1) / cv(i, 0) + 0.25 * (du(i, 1) + du(i + 2, 1)) * CuspLimiter(du(i, 1), du(i + 2, 1));
        ls(i, 2) = p(i) + 0.25 * (du(i, 2) + du(i + 2, 2)) * CuspLimiter(du(i, 2), du(i + 2, 2));

        rs(i, 0) = cv(i + 1, 0) / area(i + 1) - 0.25 * (du(i, 0) + du(i + 2, 0)) * CuspLimiter(du(i, 0), du(i + 2, 0));
        rs(i, 1) = cv(i + 1, 1) / cv(i + 1, 0) - 0.25 * (du(i, 1) + du(i + 2, 1)) * CuspLimiter(du(i, 1), du(i + 2, 1));
        rs(i, 2) = p(i + 1) - 0.25 * (du(i, 2) + du(i + 2, 2)) * CuspLimiter(du(i, 2), du(i + 2, 2));
    }
}

void flow::FluxCusp(const vec1d<double> &area) noexcept
{
}

void flow::LRStateRoe(const vec1d<double> &area, const vec1d<double> &vol) noexcept
{
    using inp::iorder;
    using inp::limfac;
    // using inp::p_ref;
    // using inp::rho_ref;
    // using inp::u_ref;
    // using inp::vol_ref;
    vec2d<double> &du = buffer2d;

    // 限制器中的无量纲所要求的
    double limfac3 = limfac * limfac * limfac;
    double rvol_ref = 1.0 / pow(vol_ref, 1.5);
    double eps2rho = limfac3 * rho_ref * rho_ref * rvol_ref;
    double eps2u = limfac3 * u_ref * u_ref * rvol_ref;
    double eps2p = limfac3 * p_ref * p_ref * rvol_ref;
    double deltr[3], deltl[3];

    if (2 <= iorder)
    {
        // 注意错了一位 第一个代表第0个
        for (int i = 1; i < imax; ++i)
        {
            du(i, 0) = cv(i, 0) / area(i) - cv(i - 1, 0) / area(i - 1);
            du(i, 1) = cv(i, 1) / cv(i, 0) - cv(i - 1, 1) / cv(i - 1, 0);
            du(i, 2) = p(i) - p(i - 1);
        }
        du(0, 0) = du(1, 0);
        du(0, 1) = du(1, 1);
        du(0, 2) = du(1, 2);
        du(imax, 0) = du(ib2, 0);
        du(imax, 1) = du(ib2, 1);
        du(imax, 2) = du(ib2, 2);
    }
    // 4.117
    auto MUSCL3 = [](double af, double bf, double eps) {
        return (bf * (2.0 * af * af + eps) + af * (bf * bf + 2.0 * eps)) /
               (2.0 * af * af + 2.0 * bf * bf - af * bf + 3.0 * eps + 1.0e-30);
    };

    double volt{0.0};
    if (iorder == 3)
    {
        for (int i = 0; i < ib2; ++i)
        {
            volt = pow(0.5 * (vol(i + 1) + vol(i)), 1.5);
            deltr[0] = MUSCL3(du(i + 2, 0), du(i + 1, 0), eps2rho * volt);
            deltl[0] = MUSCL3(du(i + 1, 0), du(i, 0), eps2rho * volt);

            deltr[1] = MUSCL3(du(i + 2, 1), du(i + 1, 1), eps2u * volt);
            deltl[1] = MUSCL3(du(i + 1, 1), du(i, 1), eps2u * volt);

            deltr[2] = MUSCL3(du(i + 2, 2), du(i + 1, 2), eps2p * volt);
            deltl[2] = MUSCL3(du(i + 1, 2), du(i, 2), eps2p * volt);

            rs(i, 0) = cv(i + 1, 0) / area(i + 1) - 0.5 * deltr[0];
            rs(i, 1) = cv(i + 1, 1) / cv(i + 1, 0) - 0.5 * deltr[1];
            rs(i, 2) = p(i + 1) - 0.5 * deltr[2];

            ls(i, 0) = cv(i, 0) / area(i) + 0.5 * deltl[0];
            ls(i, 1) = cv(i, 1) / cv(i, 0) + 0.5 * deltl[1];
            ls(i, 2) = p(i) + 0.5 * deltl[2];
        }
    }
    else
    {
        std::cerr << "roe scheme now only suport 3order\n";
        exit(EXIT_FAILURE);
    }
}

void flow::FluxRoe(const vec1d<double>& area) noexcept
{
    using inp::gam1;
    using inp::gamma;
    using inp::epsentr;
    double avs{0.0};
    double avflux[3], lrho, lu, lp, lq, lh, rrho, ru, rp, rq, rh;
    double rhot, ut, ht, ct, q2t, c2t, tmp, fdiss[3];
    double lambda[3], delta, dv, eps1, eps2, eps3;
    vec2d<double>& flux = buffer2d;

    auto EntropyCorrection = [](double lam, double dlt) {
        return (dlt < lam) ? lam : (lam*lam + dlt*dlt)/dlt*0.5;
    };

    for(int i = 0; i < ib2; ++i)
    {
        // 面积为单元中心的面积，需要求个平均
        avs = 0.5*(area(i) + area(i+1));

        lrho = ls(i, 0);
        lu = ls(i, 1);
        lp = ls(i, 2);
        lq = lu*lrho;
        // 总颔
        lh = gamma/gam1*lp/lrho + 0.5*lu*lu;

        rrho = rs(i, 0);
        ru =   rs(i, 1);
        rp =   rs(i, 2);
        rq = ru*rrho;
        rh = gamma/gam1*rp/rrho + 0.5*ru*ru;
        // average of left and right convective fluxes
        avflux[0] = lq + rq;
        avflux[1] = lq*lu + lp + rq*ru + rp;
        avflux[2] = lq*lh + rq*rh;

        // dissipative fluxes
        rhot = sqrt(lrho*rrho);
        tmp = rhot/lrho;
        ut = (lu + ru*tmp)/(1.0 + tmp);
        ht = (lh + rh*tmp)/(1.0 + tmp);
        q2t = ut*ut;
        c2t = gam1*(ht - 0.5*q2t);
        ct = sqrt(c2t);

        lambda[0] = std::abs(ut - ct);
        lambda[1] = std::abs(ut);
        lambda[2] = std::abs(ut + ct);
        delta = epsentr*ct;
        lambda[0] = EntropyCorrection(lambda[0], delta);
        lambda[1] = EntropyCorrection(lambda[1], delta);
        lambda[2] = EntropyCorrection(lambda[2], delta);
        dv = ru - lu;
        eps1 = lambda[0]*(rp -lp - rhot*ct*dv)/c2t*0.5;
        eps2 = lambda[1]*(rrho - lrho - (rp - lp)/c2t);
        eps3 = lambda[2]*(rp-lp + rhot*ct*dv)/c2t*0.5;

        fdiss[0] = eps1 + eps2 + eps3;
        fdiss[1] = eps1*(ut - ct) + eps2*ut + eps3*(ut+ct);
        fdiss[2] = eps1*(ht - ct*ut) + eps2*q2t*0.5 + eps3*(ht + ct*ut);
        // 面总的flux
        flux(i, 0) = 0.5*avs*(avflux[0] - fdiss[0]);
        flux(i, 1) = 0.5*avs*(avflux[1] - fdiss[1]);
        flux(i, 2) = 0.5*avs*(avflux[2] - fdiss[2]);
    }

    // 收集通量得到残差
    for(int i = 1; i < ib2; ++i)
    {
        rhs(i, 0) = flux(i, 0) - flux(i-1, 0);
        rhs(i, 1) = flux(i, 1) - flux(i-1, 1);
        rhs(i, 2) = flux(i, 2) - flux(i-1, 2);
    }
}